#!/usr/bin/perl
# 2011/02/06 Merge Significant peaks

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 6 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [score file] -t [threshold] -o [output file]\n";
}

my %opt;
getopts("i:t:o:", \%opt);
my $FILE_score = $opt{i};
my $FILE_output = $opt{o};
my $threshold = $opt{t};

#---------------------------------------
# threshold以上の領域を見つける
#---------------------------------------
my $fh_in;
if($FILE_score =~ /\.gz/){
	$fh_in = IO::File->new("gzip -dc $FILE_score |") or die "cannot open $FILE_score: $!";
}else{
	$fh_in = IO::File->new($FILE_score) or die "cannot open $FILE_score: $!";
}
my %tag_treatment;
my %tag_control;
my @regions;
my $pre_chr = 'I';
my $start = 0;
my $end = 0;
my $Flag_in = 0;
my $MAX_treatment = 0;
my $MAX_control = 0;
my $MAX_position = 0;
while($_ = $fh_in->getline()){
	s/\r?\n//;
	if(m/^#/){
		next;
	}

	my ($chr, $pos, $treatment, $control, $substract_nor, $logPval, $fold_change) = split /\t/;

	# tagの情報を取得
	$tag_treatment{$chr}{$pos} = $treatment;
	$tag_control{$chr}{$pos} = $control;

	# chromosomeが変わったらすべてリセット
	if($pre_chr ne $chr){
		if(($end - $start) > 99){
			push @regions, [$pre_chr, $start, $end, $MAX_position, $MAX_treatment, $MAX_control];
		}
		$start = 0;
		$end = 0;
		$MAX_position = 0;
		$MAX_treatment = 0;
		$MAX_control = 0;
	}

	if($treatment > $threshold){
		if($Flag_in == 0){
			$start = $pos;
			$Flag_in = 1;
			$MAX_position = $pos;
			$MAX_treatment = $treatment;
			$MAX_control = $control;
		}else{
			if($MAX_treatment < $treatment){
				$MAX_position = $pos;
				$MAX_treatment = $treatment;
				$MAX_control = $control;
			}
		}
		$end = $pos;
	}else{
		if(($end - $start) > 99){
			push @regions, [$pre_chr, $start, $end, $MAX_position, $MAX_treatment, $MAX_control];
		}
		$start = 0;
		$end = 0;
		$MAX_position = 0;
		$MAX_treatment = 0;
		$MAX_control = 0;
		$Flag_in = 0;
	}
	$pre_chr = $chr;
}

# すべて読み終わったら最後に溜まっている部分を出力
if(($end - $start) > 99){
	push @regions, [$pre_chr, $start, $end, $MAX_position, $MAX_treatment, $MAX_control];
}
$fh_in->close();


#---------------------------------------
# もし２つの領域が200bp以内にあったら結合させる
#---------------------------------------
my @regions2;
my $pre_start;
my $pre_end;
my $pre_MAX_pos;
my $pre_treat;
my $pre_control;
($pre_chr, $pre_start, $pre_end, $pre_MAX_pos, $pre_treat, $pre_control) = @{$regions[0]};
for(my $i = 1; $i < @regions; $i++){
	my ($chr, $start, $end, $MAX_pos, $treatment, $control) = @{$regions[$i]};

	if($chr ne $pre_chr){
		push @regions2, [$pre_chr, $pre_start, $pre_end, $pre_MAX_pos, $pre_treat, $pre_control];
		$pre_chr = $chr;
		$pre_start = $start;
		$pre_end = $end;
		$pre_MAX_pos = $MAX_pos;
		$pre_treat = $treatment;
		$pre_control = $control;
	}else{
		if(($start - $pre_end) < 200){
			$pre_end = $end;
			if($pre_treat < $treatment){
				$pre_MAX_pos = $MAX_pos;
				$pre_treat = $treatment;
				$pre_control = $control;
			}
		}else{
			push @regions2, [$pre_chr, $pre_start, $pre_end, $pre_MAX_pos, $pre_treat, $pre_control];
			$pre_start = $start;
			$pre_end = $end;
			$pre_MAX_pos = $MAX_pos;
			$pre_treat = $treatment;
			$pre_control = $control;
		}
	}
}
push @regions2, [$pre_chr, $pre_start, $pre_end, $pre_MAX_pos, $pre_treat, $pre_control];
undef @regions;



#---------------------------------------
# P-valueの計算用に候補領域の合計tagを求めて、平均値を出力
#---------------------------------------
my $fh_out = IO::File->new($FILE_output, 'w') or die "cannot write $FILE_output: $!";
$fh_out->print("#chr\tstart\tend\tmax_pos\ttreatment\tcontrol\tAVEreadPerBp__treatment\treadPerBp_control\n");
for(my $i = 0; $i < @regions2; $i++){
	my ($chr, $start, $end, $MAX_pos, $treatment, $control) = @{$regions2[$i]};

	my $SUM_treatment = 0;
	my $SUM_control = 0;
	for(my $k = $start; $k <= $end; $k++){
		$SUM_treatment += $tag_treatment{$chr}{$k};
		$SUM_control += $tag_control{$chr}{$k};
	}

	# bp辺りの平均値にする
	my $length = $end - $start + 1;
	$SUM_treatment = $SUM_treatment / $length;
	$SUM_control = $SUM_control / $length;

	$fh_out->print("$chr\t$start\t$end\t$MAX_pos\t$treatment\t$control\t$SUM_treatment\t$SUM_control\n");
}







